library(readxl)
## Warning: package 'readxl' was built under R version 4.1.3
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(RColorBrewer)
#this.dir = rstudioapi::getActiveDocumentContext()$path
#setwd(dirname(this.dir))
linaclotide.data = read_xls(dir(pattern='XLS'), skip = 4,
                      col_types = c(
                        'text', 'skip', 'numeric', 'skip', 'text', 'text',
                        'numeric', 'skip', 'numeric', 'numeric'
                      ))
## Warning: Coercing text to numeric in H7 / R7C8: '2.442'
## Warning: Coercing text to numeric in J7 / R7C10: '172.760'
## Warning: Coercing text to numeric in K7 / R7C11: '345.369'
## Warning: Coercing text to numeric in H8 / R8C8: '2.438'
## Warning: Coercing text to numeric in J8 / R8C10: '172.931'
## Warning: Coercing text to numeric in K8 / R8C11: '345.614'
## Warning: Coercing text to numeric in H9 / R9C8: '2.440'
## Warning: Coercing text to numeric in J9 / R9C10: '172.148'
## Warning: Coercing text to numeric in K9 / R9C11: '345.825'
## Warning: Coercing text to numeric in H10 / R10C8: '2.443'
## Warning: Coercing text to numeric in J10 / R10C10: '0.602'
## Warning: Coercing text to numeric in K10 / R10C11: '1.512'
## Warning: Coercing text to numeric in H11 / R11C8: '2.442'
## Warning: Coercing text to numeric in J11 / R11C10: '110.138'
## Warning: Coercing text to numeric in K11 / R11C11: '213.756'
## Warning: Coercing text to numeric in H12 / R12C8: '2.446'
## Warning: Coercing text to numeric in J12 / R12C10: '54.947'
## Warning: Coercing text to numeric in K12 / R12C11: '102.733'
## Warning: Coercing text to numeric in H13 / R13C8: '2.450'
## Warning: Coercing text to numeric in J13 / R13C10: '12.208'
## Warning: Coercing text to numeric in K13 / R13C11: '21.988'
## Warning: Coercing text to numeric in H14 / R14C8: '2.450'
## Warning: Coercing text to numeric in J14 / R14C10: '6.031'
## Warning: Coercing text to numeric in K14 / R14C11: '10.765'
## Warning: Coercing text to numeric in H15 / R15C8: '2.451'
## Warning: Coercing text to numeric in J15 / R15C10: '1.347'
## Warning: Coercing text to numeric in K15 / R15C11: '2.403'
## Warning: Coercing text to numeric in H16 / R16C8: '2.452'
## Warning: Coercing text to numeric in J16 / R16C10: '0.660'
## Warning: Coercing text to numeric in K16 / R16C11: '1.173'
## Warning: Coercing text to numeric in H17 / R17C8: '2.453'
## Warning: Coercing text to numeric in J17 / R17C10: '0.154'
## Warning: Coercing text to numeric in K17 / R17C11: '0.267'
## Warning: Coercing text to numeric in H18 / R18C8: '2.445'
## Warning: Coercing text to numeric in J18 / R18C10: '109.991'
## Warning: Coercing text to numeric in K18 / R18C11: '213.800'
## Warning: Coercing text to numeric in H19 / R19C8: '2.447'
## Warning: Coercing text to numeric in J19 / R19C10: '57.127'
## Warning: Coercing text to numeric in K19 / R19C11: '106.836'
## Warning: Coercing text to numeric in H20 / R20C8: '2.452'
## Warning: Coercing text to numeric in J20 / R20C10: '12.077'
## Warning: Coercing text to numeric in K20 / R20C11: '21.694'
## Warning: Coercing text to numeric in H21 / R21C8: '2.452'
## Warning: Coercing text to numeric in J21 / R21C10: '6.009'
## Warning: Coercing text to numeric in K21 / R21C11: '10.728'
## Warning: Coercing text to numeric in H22 / R22C8: '2.453'
## Warning: Coercing text to numeric in J22 / R22C10: '1.245'
## Warning: Coercing text to numeric in K22 / R22C11: '2.233'
## Warning: Coercing text to numeric in H23 / R23C8: '2.455'
## Warning: Coercing text to numeric in J23 / R23C10: '0.623'
## Warning: Coercing text to numeric in K23 / R23C11: '1.108'
## Warning: Coercing text to numeric in H24 / R24C8: '2.453'
## Warning: Coercing text to numeric in J24 / R24C10: '0.135'
## Warning: Coercing text to numeric in K24 / R24C11: '0.216'
## Warning: Coercing text to numeric in H25 / R25C8: '2.446'
## Warning: Coercing text to numeric in J25 / R25C10: '101.965'
## Warning: Coercing text to numeric in K25 / R25C11: '195.745'
## Warning: Coercing text to numeric in H26 / R26C8: '2.447'
## Warning: Coercing text to numeric in J26 / R26C10: '53.896'
## Warning: Coercing text to numeric in K26 / R26C11: '100.139'
## Warning: Coercing text to numeric in H27 / R27C8: '2.451'
## Warning: Coercing text to numeric in J27 / R27C10: '11.378'
## Warning: Coercing text to numeric in K27 / R27C11: '20.314'
## Warning: Coercing text to numeric in H28 / R28C8: '2.453'
## Warning: Coercing text to numeric in J28 / R28C10: '5.725'
## Warning: Coercing text to numeric in K28 / R28C11: '10.170'
## Warning: Coercing text to numeric in H29 / R29C8: '2.453'
## Warning: Coercing text to numeric in J29 / R29C10: '1.226'
## Warning: Coercing text to numeric in K29 / R29C11: '2.200'
## Warning: Coercing text to numeric in H30 / R30C8: '2.453'
## Warning: Coercing text to numeric in J30 / R30C10: '0.589'
## Warning: Coercing text to numeric in K30 / R30C11: '1.043'
## Warning: Coercing text to numeric in H31 / R31C8: '2.453'
## Warning: Coercing text to numeric in J31 / R31C10: '0.142'
## Warning: Coercing text to numeric in K31 / R31C11: '0.246'
## Warning: Coercing text to numeric in H32 / R32C8: '2.446'
## Warning: Coercing text to numeric in J32 / R32C10: '101.758'
## Warning: Coercing text to numeric in K32 / R32C11: '194.768'
## Warning: Coercing text to numeric in H33 / R33C8: '2.446'
## Warning: Coercing text to numeric in J33 / R33C10: '107.074'
## Warning: Coercing text to numeric in K33 / R33C11: '206.443'
## Warning: Coercing text to numeric in H34 / R34C8: '2.445'
## Warning: Coercing text to numeric in J34 / R34C10: '101.591'
## Warning: Coercing text to numeric in K34 / R34C11: '194.717'
## Warning: Coercing text to numeric in H35 / R35C8: '2.452'
## Warning: Coercing text to numeric in J35 / R35C10: '11.976'
## Warning: Coercing text to numeric in K35 / R35C11: '21.640'
## Warning: Coercing text to numeric in H36 / R36C8: '2.451'
## Warning: Coercing text to numeric in J36 / R36C10: '12.022'
## Warning: Coercing text to numeric in K36 / R36C11: '21.398'
## Warning: Coercing text to numeric in H37 / R37C8: '2.451'
## Warning: Coercing text to numeric in J37 / R37C10: '11.688'
## Warning: Coercing text to numeric in K37 / R37C11: '20.788'
## Warning: Coercing text to numeric in H38 / R38C8: '2.452'
## Warning: Coercing text to numeric in J38 / R38C10: '0.462'
## Warning: Coercing text to numeric in K38 / R38C11: '0.872'
## Warning: Coercing text to numeric in H39 / R39C8: '2.453'
## Warning: Coercing text to numeric in J39 / R39C10: '0.398'
## Warning: Coercing text to numeric in K39 / R39C11: '0.711'
## Warning: Coercing text to numeric in H40 / R40C8: '2.454'
## Warning: Coercing text to numeric in J40 / R40C10: '0.388'
## Warning: Coercing text to numeric in K40 / R40C11: '0.709'
## Warning: Coercing text to numeric in H41 / R41C8: '2.455'
## Warning: Coercing text to numeric in J41 / R41C10: '0.119'
## Warning: Coercing text to numeric in K41 / R41C11: '0.220'
## Warning: Coercing text to numeric in H42 / R42C8: '2.454'
## Warning: Coercing text to numeric in J42 / R42C10: '0.120'
## Warning: Coercing text to numeric in K42 / R42C11: '0.212'
setClass('sample.prep',
         slots = list(compound='character', n.cal='numeric', reps='numeric'))
linaclotide.prep = new('sample.prep', compound='linaclotide', n.cal=7, reps=3)

n = linaclotide.prep@n.cal
lin.data.1 = linaclotide.data %>%
    slice(-c(1:4)) %>%
    cbind('stock' = c(rep(1:3, each = 7), rep(seq(4, 6), 4), 'blank')) %>%
    cbind('type' =
              c(rep('Cal', n * linaclotide.prep@reps),
                rep('QC', 12),
                'blank')) %>%
    cbind('concentration' = c(rep(c(100, 50, 10, 5, 1, 0.5, 0.1), 3),
                              rep(c(100, 10, 0.3, 0.1), each = 3),
                              0)) %>%
    rename(Area = `Peak Area`) %>%
    replace(is.na(.), 0)

lin.data.cal = lin.data.1 %>% filter(type == 'Cal')
cal.curve = lm(Area ~ concentration, data = lin.data.cal, weights = 1/concentration)
summary(cal.curve)
## 
## Call:
## lm(formula = Area ~ concentration, data = lin.data.cal, weights = 1/concentration)
## 
## Weighted Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.20650 -0.12310  0.06483  0.26783  0.59900 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.05765    0.06429   0.897    0.381    
## concentration  2.07752    0.01818 114.249   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.402 on 19 degrees of freedom
## Multiple R-squared:  0.9985, Adjusted R-squared:  0.9985 
## F-statistic: 1.305e+04 on 1 and 19 DF,  p-value: < 2.2e-16
cal.plot.lin = ggplot(lin.data.1, mapping = aes(concentration, Area)) +
    geom_point(aes(colour = stock, shape = type)) +
#    stat
    geom_smooth(data = lin.data.cal, 
        method = 'lm', mapping = aes(weight = 1/concentration), 
        formula = y ~ x, ) +
    labs(title='linaclotide calibration', x='concentration, ug/mL',
         y='area, mAU*s') +
    theme_classic()
#cal.plot
#
#cal.plotly = plot_ly(data = data.group, x = ~concentration, y = ~Area,
#                     color = ~stock, colors = "Reds", symbol = ~type)
#cal.plotly
cal.plot.lin

ggplotly(cal.plot.lin)
data.post.lin = lin.data.1 %>%
    group_by(stock) %>%
    mutate(calc.conc = (Area - cal.curve$coefficients[1])/
               cal.curve$coefficients[2],
            accuracy = calc.conc / concentration)
data.post.lin
data.post.lin.stock = data.post.lin %>% summarise(avg.acc = mean(accuracy))
data.post.lin.stock